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ABSTRACT 
A  semi-implicit  method  for  solving  the  three-dimensional 
magnetohydrodynamic  equations  on  long  time  scales  is  presented. 
Standard  explicit  methods  must  use  time  steps  which  are  constrained  by 
a  Courant-Friedrichs-Lewy  condition  due  to  the  fast  compressional  and 
shear  Alfven  motion.  This  semi-implicit  method  eliminates  both  of 
these  restrictions  so  that  very  large  time  steps  are  permitted.  The 
method  is  simple  to  implement  and  the  computation  time  for  one  time 
step  is  essentially  the  same  as  for  explicit  methods.  Numerical  test 
results   in  slab  and   cylindrical   geometry  are  presented. 
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I.   INTRODUCTION 

The  magnetohydrodynamic  (MHD)  equations  are  used  extensively  to 

1  2 
Study   the   macroscopic  behavior  of  plasmas.  •   Three-dimensional 

time-dependent  computations  are  difficult  due  to  the  presence  of  widely 

disparate  time  scales   in  the  MHD  equations.   Explicit  methods  are 

forced  to   use   time   steps   limited   by   a   very   restrictive 

Courant-Friedrichs-Lewy   (CFL)   condition   imposed   by   the   fast 

compressional  (magnetosonic)  motion.   Implicit   schemes   for   the 

three-dimensional  MHD  equations  have  been  developed^  but  have  generally 

had  the  disadvantages  of  being  very  complicated  to  implement  and 

requiring  the  solution  of  large  block  matrix  equations.   A  common 

approach  to  eliminate  the  fast  compressional  CFL  condition  is  to  make 

analytic  approximations  such  as  in  the  reduced  equations    (inverse 

aspect  ratio  expansion)  or  incompressible  models.   Unfortunately,   in 

many   problems   (e.g.    incompressible  reversed-f ield  pinch  dynamo 

computations'^)  important  features  of  the  physical   system  are  also 

eliminated  by  these  approximations.   Recently,  two  different  methods 

have  been  developed   in  order  to  solve  the  full   compressible  MHD 

equations,   but  without  a  CFL  time  step  restriction  due  to  the  fast 

modes.   One  approach  is  that  of  Aydemir  and  Barnes, °  where  an   implicit 

pressure  advance   is  used  to  eliminate  the  fast  mode  constraint.   The 

QIC 
other  method  is  the  semi-implicit  method  of  Harned  and  Kerner^'    in 

which  new  terms  are  added  to  the  time  discretized  MHD  equations.   These 

new  terms  do  not  affect  the  solution  as  At  -»  0,   but   still   produce  a 

method  that  is  unconditionally  stable  with  respect  to  the  fast  modes. 

Once   the  fast  compressional  time  step  constraint  has  been 

eliminated,  the  time  step   is  then  limited  by  a  shear  Alfven  CFL 
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constraint.  This  constraint  is  particularly  severe  in  the  case  of  the 
reversed-f ield  pinch,  in  which  the  toroidal  and  poloidal  magnetic 
fields  are  comparable.  In  tokamak  plasmas  the  shear  Alfven  CFL 
condition  becomes  very  restrictive  in  the  nonlinear  phase  of  resistive 
instabilities.  The  methods  of  Refs.  8-10,  as  well  as  incompressible 
and  reduced  equation  methods,  all  require  a  substantial  reduction  in 
time  step  in  the  highly  nonlinear  stages  of  resistive  instabilities, 

which  ultimately  limit  their  capabilities.   In  this  paper  we  generalize 

Q  1  n 
the  method  of  Harned  and  Kerner^'   to  remove  the  shear  Alfven  CFL  time 

step  restriction. 

The  method  to  be  described  here  is  designed  to  eliminate  the  shear 

Alfven  time  step  restriction  as  well  as  that  due  to  the  fast 

compressional  waves.   It  is  much  simpler  to  implement  than  an  implicit 

method  and  the  computation  time  required  for  a  time  step  is  essentially 

the  same  as  that  for  an  explicit  method.   Section  II  describes  the 

model  and  the  semi-implicit  method.   The  results  of  numerical  tests  are 

presented  in  Sec.   Ill  and  conclusions  are  given  in  Sec.   IV. 


II.   SEMI-IMPLICIT  MHD  METHOD 

The  compressible  resistive  MHD  equations  may  be  written 


ll  -  -  (v.V)v  +  1  [(VxB)  X  B-VP]  (1) 

3t  p 


H  -  VxB  -  nVxVxA  (2) 

at 


|£  =  -  V.(pv)  (3) 

at 


3P     ~  ~ 

I_  =  -  vVP  -  YPV'V  +  dissipative  terms  (M) 

9t 

B  =  VxA  (5) 


where  v  is  the  velocity,  B  the  magnetic  field,  P  the  plasma  pressure,  p 
the  density  and  n  the  resistivity.  A  is  the  vector  potential  and  the 
gauge  condition  4)  =  0  has  been  used.  Resistive  instabilities  are 
important  in  the  analysis  of  fusion  plasmas.  However,  they  evolve  on  a 
time  scale  which  is  long  compared  to  ideal  time  scales  (e.g.,  fast 
compressional  and  shear  Alfven).  Therefore,  in  order  to  solve  Eqs. 
(1)-(5)  numerically,  it  is  desirable  to  develop  a  method  that  is  not 
forced  to  use  time  steps  constrained  by  the  rapid  ideal  motion. 

To  remove  the  tiraestep  constraint  imposed  by  the  fast 
compressional  modes  a  semi-implicit  method  was  used  in  Refs.  9  and  10. 
Because  of  the  flexibility  of  this  method  it  is  natural  to  try  to 
extend  it  to  eliminate  the  shear  Alfven  time  step  constraint  as  well. 
In  a  semi-implicit  method  one  advances  all  of  the  terms  in  the  original 
equations  explcitly.   New  terms  are  added  to  the  time  discretized 
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equations  which  do  not  affect  the  solution  in  the  limit  At  ->•  0   (i.e. 
the  method   is   still  consistent).   Then  some  of  these  "semi-implicit" 
terms  are  treated   implicitly.    As  an  example   consider   the  simple 
hyperbolic  system, 


8u     9v  ,., 

=  a  -—  (6) 


3t     9x 


9v     3u  ,„. 

=  a  -_  .  (7) 


3t     9x 


These  equations  may  be  rewritten 


il-   a^i^   .  (8) 

at^    3x^ 


To  develop  a  semi-implicit  method,  a  new  term  is  subtracted  from  each 
side  of  Eq.   (8), 


3t2      3x^      3x2      3^2 


a^  is  a  constant  coefficient  and  will  be  chosen  primarily  from 
numerical  stability  considerations.  Equation  (9)  is  then  time 
differenced  as 


u"^!  -  a2(At)2  ifl^  =  u"  .  At(i^]"  *  a2(At)2  l!^  -  a2(At)2  ^      (10) 
°       3x2  3t  9,2     o      9,2 


A  practical  way  to  time  advance  Eqs.   (6)  and  (7)  in  order  to  obtain  an 
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algorithm    equivalent   to   Eq.    (10)   is   with   a   two   step 
predictor- corrector 


u*  -  u"  *  caAt  ^  (11) 

8x 


V*  -  v"  -  aaAt  ^  (12) 

dx 

un*1  -  a2(At)2  i^   =  u"  wAt  ^  -  a2(At)2  l!l^       (13) 
8x2  9x  3^2 


v"*1  =  v"  *  (a/2)At  [S}i.   +  i!^ ]  (U) 

3x     3x 


with  0.5  ^  a  ^  1.0.  This  method  is  unconditionally  stable  if 
a^  >  (a^/l 6)( 1 +2a)2.  For  a  simple  one-dimensional  case  like  this  the 
semi-implicit  method  offers  no  advantages  over  an  implicit  method. 
However,  if  a  problem  is  three-dimensional  and  spectral  in  two 
dimensions  with  nonconstant  coefficients,  an  implicit  method  is 
difficult  because  it  requires  performing  convolutions  in  the  implicit 
terms.  This  leads  to  a  complicated  matrix  equation.  In  the 
semi-implicit  method  the  constant  coefficient  semi-implicit  terms 
require  no  convolutions;  therefore,  only  a  simple  tridiagonal  solution 
is  needed.  Another  important  feature  of  the  method  is  the  use  of  the 
second  order  equation  for  discretization.  If  the  method  was  applied  to 
the  original  first  order  equations,  the  sign  of  a  becomes  important. 
If  a  changes  sign,  the  constant  coefficient  semi-implicit  term  would 
not  stabilize  the  system  because  in  the  region  where  a  and  a^  are  of 
opposite  sign  the  semi-implicit  term  would  actually  be  destabilizing. 
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However,   in  the  second  order  equation,  Eq.   (9),  the  coefficient  only 
enters  as  a^  so  that  the  sign  does  not  matter  and  the  semi-implicit 
treatment  becomes  very  simple. 

In  Refs.  9  and  10  the  fast  compressional  time  step  constraint  was 
eliminated  by  choosing  a  semi-implicit  term  having  the  same  form  as  the 
fast  modes.  This  term  was  subtracted  from  each  side  of  a  second  order 
equation.  Then  the  MHD  equations,  with  the  new  semi-implicit  terms 
included,  were  differenced  in  a  predictor-corrector  scheme  analogous  to 
Eqs.  (11)-(14).  To  extend  the  method  to  also  eliminate  the  shear 
Alfven  time  scale,  we  must  derive  a  new  semi-implicit  term.  First,  the 
MHD  equations  are  linearized,  assuming  a  uniform  magnetic  field, 
density,  and  pressure.  Then,  after  some  algebra,  Eqs.  (l)-(5)  reduce 
to 


a^v      -  -   - 

L-L   =  [VxVx(vxBq)]xBq  +  YPqV(V.v)  .  (15) 

3t2 


To  arrive  at  a  semi-implicit  term,   B^  is  replaced  with  a  vector 

—  y^  /S  ^ 

quantity  with   constant   coefficients,  C^  =  C  x  =  C  y  +  C^z,  where  C^, 

C  ,  and  C  are  all  constant  in  space.   The  pressure   term   is  dropped, 
y     2 

because  it  enters  in  the  magnetosonic  modes  with  the  same  form  as  the 
perpendicular  magnetic  field.  Therefore,  the  semi-implicit  term  is 
just 


[VxVx(vxCq)]  X  Cq  . 


One  now  naively  hopes  that  this  term  could  be  used  in  a  predictor- 


-8- 

corrector  method  to  eliminate  the  shear  Alfven  time  step  constraint   in 
the  same  way  as   the   fast  mode  time  step  constraint  was  eliminated. 

This  philosophy  suggests  the  following  algorithm.  First,   a  simple 
explicit  predictor  advance  is  performed, 

V*  =  v"  +  aAt  -1  F  (p.v.  B,  P)"  (16) 
p" 

A*  -  a"  *  oAt  v"xB"  (17) 

P*  =  P"-  aAt(v"-VP"  +  YP"V.v")  (18) 

p*  »  p"-  cxAtV.(pv)"  (19) 

B*  =  VxA*  (20) 

The  semi-implicit  term  is  included  in  the  corrector  velocity  advance: 


^n*1  _  UUf  [v.Vx(v"*^Co)]xCo  =  v"  +  ^  F(p.v,B,P)*  -  lAl^  [VxVx(>xCq)  ]  .C^   (2 
p*  p  p 


The  new  velocity  is  used  in  the  magnetic  field,  pressure,  and  density 
advances, 


?*  =  A,  -  i^  C(v""Uv")xB*]  (22) 

n    2 

pn+1  =  pM  _  At  ^(;n+1^;n).^p«  ^  yPV- ( v"^^ +>)]  (23) 

2 


pn-1  =  pn  _  At  v.[p*(v"*^v")]  (2^) 
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Finally,  a  time  step  is  completed  with  a  semi-implicit  resistive 
advance 

A"*^  +  AtriQVxVxA"*^  •=  A**  -  AtnVxVxA**  +  AtngV^VxA**  (25) 

The  resistive  part  of  the  MHD  equations  has  been  split  from  the  ideal 
part  and  treated  semi-implicitly  so  that  the  resistive  advance  does  not 
impose  an  additional  timestep  constraint.  As  in  the  Harned-Kerner 
algorithm J"^  the  treatment  is  semi-implicit  rather  than  implicit  so 
that  nonuniform  resistivity  may  be  used  without  requiring  convolutions 
to  be  performed  in  the  implicit  part. 

It  is  hoped  that  if  a  set  of  conditions  like  C  >  B^,  C  >  B  and 
C  >  B  is  satisfied,  the  method  would  be  unconditionally  stable  with 
respect  to  both  the  fast  compressional  and  shear  Alfven  modes.  For  a 
one-dimensional  case  a  linear  stability  analysis  shows  that  this  is 
true.  However,  the  critical  difference  between  this  method  and  that 
used  for  the  fast  modes  is  that  C  is  a  vector  rather  than  a  scalar. 
Unfortunately,  one  finds  for  a  two-dimensional  stability  analysis  that 
unconditional  stability  for  the  fast  and  shear  modes  can  occur  only 
when  C  is  parallel  to  B.  This  means  that  when  B  is  nonuniform,  C^ 
must  be  nonuniform  as  well.  This  makes  the  method  equivalent  to  an 
implicit  method,  requiring  the  solution  of  complicated  matrix  equations 
for  the  time  advance.  Such  a  method  is  impractical  and  defeats  the 
purpose  of  the  semi-implicit  method. 

To  solve  this  problem,  it  is  important  to  realize  that  the  form  of 
the  semi-implicit  terms  is  completely  arbitrary  (although  extreme 
choices  may  seriously  degrade  the  accuracy  of  the  method).    The   terms 
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3hould  be  chosen  to  enhance  the  numerical  stability  of  the  method,  yet 
still  be  simple  to  treat  implicitly.  In  a  two-dimensional  stability 
analysis,  the  terms  having  semi-implicit  coefficients  C^Cj,  with  i  »<  J, 
are  destabilizing  when  C^  and  B^  are  not  parallel.  When  these  terms 
are  set  to  zero,  an  unconditionally  stable  algorithm  results.  In 
addition,  this  new  modified  algorithm  has  the  desirable  property  that 
the  semi-implicit  coefficients  enter  only  as  C^,  so  that  the  sign  of 
the  coefficients  does  not  matter. 

To  demonstrate  the  linear  stability  properties  of  the  method, 
consider  a  two-dimensional  cold  ideal  case  with  a  uniform  density  and 
equilibrium  magnetic  field.   The  MHD  equations  reduce  to 


—   =  (7xB)xB  (26) 

3t 


M  =  vxB  (27) 

3t 


Let  A  -  A  z,  V  =  v^x  ♦  v  y,  and  3/3z  =  0.0.   Then  a  method  similar  to 
that  of  Eqs.   (l6)-(24)  is 

v'  =  v"  *  At(VxB")xB"  (28) 

A*  =  A^  *  At(v"xB")2  (29) 

V"*^  -  (At)2[VxVx(v"*^CQ)]xCo  =  V"  *  At(VxB*)xB*  -  (At)2[VxVx( v"xCo)]xCq   (30) 

j^n*1  ,  ^n  ,   (vn*1xB*),  (SD 
z     z          z 

except  that  Eq.   (30)  is  modified  by  setting  the  terms  with  C^C^   equal 
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to  zero.   The  spatial  differencing  Is  computed  with  a  spectral 
representation, 


f(x)  -   I   fnme'"'^*'"^ 


m.n 


The  system,    Eqs.      (28)-(31  )    with  the  C^^C     terms   dropped,    may  be  'written 
as  the  matrix  equation 


1*n2c2 


0 

AtB, 


AtB^ 


0 

—                              Si 

^x 

n+1 

0 

^ 

- 

1 

/z_ 

1*n2(C^-b2)         BjjByN^ 


B^B 


/' 


-ByN^/At 


1+n2(C^-B^)         BjjN^/At 


'  X 

^y 
'x 


n 


where  N^   is     defined     as     N^  =    (m^+n   )(At)    .        Equation      (32)      is     then 
rewritten   in  the  form 


(3: 


r  nn+i 


n 


(33) 


where  G  is  the  amplification  matrix, 


G  = 


1-Y 


Byd-X-Y) 


y 

1-X 


-Bj^d-X-Y) 


-Y/B, 


X/B. 


1-X-Y 


(3^) 
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In  Eq.  {3k)   X  and  Y  are  defined  by 


X  _  N^B^/d+N^C^) 


and 


Y  _  N^B^/d+N^C^) 


To  determine  the  linear  stability  of  the  algorithm,  the  eigenvalues  of 
the  amplification  matrix  are  computed.  In  addition  to  the  trivial 
result,  (jo  =  1  ,  they  are 


0)  =  1  -  Z  ±  /Z(Z-D  (35) 

where  Z  =  X  +  Y.  When  Z  <  1.0  the  eigenvalues  lie  inside  the  unit 
circle.  Therefore,  as  desired,  the  method  is  unconditionally  stable  as 
long  as  C^  >  B„  and  C„  >  B„.  This  is  not  the  case  if  the  C  C  cross 
terms  are  retained  in  the  semi-implicit  corrector.  If  smaller  values 
of  a  are  used  in  the  predictor,  the  semi-implicit  coefficients  may  also 
be  reduced. 

The  complete  semi-implicit  predictoi — corrector  algorithm  that 
eliminates  both  the  shear  Alfven  and  fast  compressional  CFL  time  step 
constraints  is  given  by  Eqs.  (l6)-(25),  with  the  crucial  modification 
that  all  terms  with  C^C.  are  replaced  by  C^C.6^^  in  Eq.  (21).  Solving 
this  system  is  not  difficult.  A  spectral  representation  is  used  in  two 
directions.  In  the  x-direction  (or  the  radial  direction  for 
cylindrical  coordinates)  centered  finite  differences  are  used.  All  of 
the  quantities  in  the  ideal  advance  are  advanced  explicitly  with  the 
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exception  of  the  velocity  corrector.  In  the  velocity  corrector  the 
semi-implicit  term  on  the  left  hand  side  is  treated  implicitly,  but  no 
convolutions  are  required.  Therefore,  this  equation  is  a  block 
tridiagonal  matrix  equation  with  3x3  blocks.  This  equation  may  be 
solved  easily  without  adding  any  significant  amount  of  computation  time 
per  time  step,  vrfien  compared  to  an  explicit  method.  This  is  because 
for  any  case  with  at  least  a  few  modes,  the  computation  time  is 
dominated  by  the  convolutions.  As  the  plasma  evolves  nonlinearly  in 
time  the  coefficients,  C^,  may  be  varied  in  time  (and  in  space,  if 
desired)  to  preserve  linear  stability  with  respect  to  the  new  values  of 
magnetic  field  and  pressure.  The  method  described  here  is  first  order 
accurate  in  time.  First  order  accuracy  is  not  an  inherent  feature  of 
semi-implicit  methods.  Second  order  accurate  schemes  simply  require 
more  complicated  semi-implicit  corrector  equations.  A  second  order 
accurate  semi-implicit  method  is  being  used  by  Joyce  for  beam 
propagation  problems  and  the  Harned-Kerner  method  of  Ref.  10  has  also 
been  implemented  in  a  second  order  accurate  version. 


III.   NUMERICAL  TESTS 

The  new  seml-impliclt  algorithm  described  in  Sec.  II  has  been 
tested  in  slab  and  cylindrical  geometry  to  demonstrate  its  numerical 
stability  properties.  In  a  slab  we  assume  a  constant  density,  p  =  1.0. 
The  y  and  z  directions  are  periodic  and  are  treated  with  a  spectral 
representation.  Finite  differences  are  used  in  the  x-direction.  Rigid 
conducting  wall  boundary  conditions  are  imposed  at  x  =  0.0  and  x  =  1.0. 
First,  a  uniform  magnetic  field  is  initialized  with  B  =  1.0  and 
B  =  0.2.  We  apply  the  following  three-dimensional  perturbation  with 
m  =  1  and  n  •=  1  ,  where  6  is  the  perturbation  amplitude: 

v^  =  (0.96/2ir)6  sin  (my  +  0.2  nz)  sin  (2ttx) 
V  =  1.0  5  cos  (my  +  0.2  nz)  cos  (2iix) 
v^  -  -0.2  6  cos  (my  +  0.2  nz)  cos  (2itx) 

This   is  an  excitation  of  a  shear  Alfven  wave   with   frequency 

0)  =  k«B  =  0.4.   We  use  41   grid  points  and  keep  20  modes  from  m  =  1, 

n  =  1  to  m  =  20,  n  =  20.   Although  the  wave  is  linear,   the  code  is 

nonlinear  so  that  all  of  the  modes  are  excited.   If  we  do  not  use  the 

semi-implicit  method,  i.e.   set  C^  =  0.0,   then  numerical   instability 

results  if  the  time  step  exceeds  the  usual  CFL  condition  on  the  fast 

modes,  At  <  Ax/v_  =  0.025.   Furthermore,  even  if  the  fast  mode  time 

step  constraint   is  eliminated,  the  shear  Alfven  time  step  constraint 

limits  the  time  step  to  At  <  0.24.   To  simulate  the  Alfven  wave  we 

first  choose  a  time  step  of  At  =  0.1.   which  is  four  times  the  normal 

CFL  limit.  The  numerical  parameter,   a,   of  the  predictor-corrector 

method   is  set  to  a  =  0.6.   We  use  C_  =  0.7  and  C^  =  0.2.   The  kinetic 

z  y 
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energy  of  the  wave  is  shown  in  Fig.   1.   The  expected  period  of  t   -  5-n 

is  produced  correctly.   In  Fig.   2  we  show  the  result  of  a  simulation 

with  At  -  0.5,  more  than  twice  the  shear  Alfven  limit.   Some  damping  of 

the  wave  kinetic  energy  is  apparent  due  to  the  predictor-corrector 

method.   The  wave  frequency  is  still  correct.   As  the  time  step  is 

increased  further,   the  method  remains  stable  but  the  wave  becomes 

poorly  represented.   In  fact,  the  time  step  can  be  raised  to  extremely 

large  values,   such  as  At  =  200,  and  the  algorithm  will  be  stable  but 

the  wave  is  immediately  damped.   We  note  that  if  the  cross  terms   (e.g. 

C  C  )   in  the  semi-implicit  terms  are  retained,  the  method  fails,  as 
X  y 

expected. 

As  a  second  test  we  use  a  uniform  magnetic  field,   B^  =  1.0,   and 
apply  a  two-dimensional  perturbation  with  . 


V  «=  6  cos  (ny)  cos  (2itx) 


Twenty  modes  with  m  =  0  are  retained  with  the  largest  being  n  =  200. 

An  additional  small  random  perturbation  is  also  applied  so  that  all  of 

the  modes  are  excited.   Setting  C^  =  0.7  and  C  =  0.0  would  be 

sufficient  to  have  an  unconditionally  stable  algorithm  for  a  linear 

wave.   However,   if  the  wave  has  a  large  amplitude,  then  B  becomes 

substantial  and  it  may  impose  a  stability  limit  as  well.    For  a  case 

when  the  perturbed  B  =  0.3  and  At  =  0.1,  the  method  goes  immediately 

numerically  unstable.   However,  when  C  =  0.3   is   used   unconditional 

stability   again  results.    The   kinetic   energy   for   the   case   with 

C  =  0.7,  C„  =  0.3,  and  At  =  0.1  is  shown  in  Fig.   3-   The  algorithm  is 
^       y 
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stable  even  with  the  large  amplitude  oscillatory  magnetic  field,  B  , 
present  and  the  correct  frequency  is  obtained. 

In  cylindrical  geometry  we  initialize  an  equilibrium  with 
B  «  1.0,  J_  -  J^(l-r^/a^),  and  k„  =  0.33.  The  safety  factor  at  the 
wall,  q(a)  -  2iTrB2(a)/LBQ(a) ,  is  set  to  q(a)  -  2.93-  The  safety  factor 
for  this  profile  is  shown  in  Fig.  4.  The  resistivity  is  set  to 
n  -  10~^.  200  radial  grid  points  are  used  and  twenty  modes  from  m  =  2, 
n  =  1  to  m  =  40,  n  =  20  are  retained.  The  time  step  is  set  to  At  =  1 . 5 
radial  Alfven  transit  times.  The  normal  CFL  limit  for  an  explicit  code 
would  be  At  <  0.005  for  the  fast  modes  and  At  <  0,12  for  the  shear 
Alfven  modes.  Hence,  both  of  these  conditions  have  been  violated  by 
more  than  a  factor  of  ten.  For  this  case  we  allow  the  semi-implicit 
coefficients  to  vary  with  radius  and  time.  At  each  time  step  the 
components  of  C  are  set  to  0.9  times  the  maximum  value  of  their 
respective  magnetic  field  components  at  a  given  radius.  A  random 
perturbation  is  applied  so  that  the  fastest  growing  unstable  mode 
should  eventually  dominate.  We  expect  to  observe  an  m  =  2,  n  =  1 
instability  with  a  growth  rate  of  Y  =  0.0046.  The  kinetic  energy  as  a 
function  of  time  is  plotted  in  Fig.  5.  The  growth  rate  obtained  from 
the  simulation  is  Y  =  0.0044,  giving  good  agreement.  The  profiles  of 
Vj,  and  J^  are  shown  in  Figs.  6  and  7,  respectively.  In  Fig.  8  are 
plotted  helical  flux  contours,  showing  the  m  =  2,  n  =  1  island  near  the 
saturation  of  linear  growth.  The  time  step  may  be  increased  further. 
However,  very  large  time  steps  degrade  the  accuracy  of  the  results.  A 
time  step  of  five  Alfven  transit  times  gives  a  substantially  smaller 
growth  rate  for  the  instability,  Y  =  0.0028,  and  the  velocity  and 
current  profiles  show  a  departure   from  the  proper  eigenfunctions. 
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Therefore,  for  this  method  the  primary  limiting  factor  on  the  time  step 
appears  to  be  accuracy  rather  than  stability.  We  also  note  for  this 
case  with  20  modes  that  the  computation  time  required  per  time  step  is 
only  5%  more  than  for  an  explicit  advance,  even  though  much  larger  time 
steps  are  permitted. 
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IV.   CONCLUSIONS 

A  new  semi-implicit  method  for  solving  the  MHD  equations  in  three 
dimensions  has  been  developed.  The  method  allows  simulations  to  use 
very  large  time  steps  since  both  the  shear  Alfven  and  the  fast 
compressional  CFL  conditions  have  been  eliminated.  The  method  is 
simple  to  implement.  The  most  difficult  part  of  the  time  advance  is 
the  solution  of  a  block  tridiagonal  system  of  3x3  blocks.  Therefore, 
the  method  requires  virtually  the  same  amount  of  computation  time  per 
time  step  as  an  explicit  method. 

The  method  has  been  tested  on  simple  problems  in  slab  and 
cylindrical  geometry,  which  verify  the  unconditional  stability  of  the 
algorithm  with  respect  to  both  shear  Alfven  and  fast  compressional 
modes.  The  primary  limiting  factor  on  time  step  now  appears  to  be 
accuracy.  In  the  algorithm  of  Ref.  10,  the  results  using  first  order 
and  second  order  accurate  methods  are  essentially  identical.  In  that 
algorithm  time  steps  were  limited  by  the  shear  Alfven  CFL  condition  and 
were  therefore  sufficiently  small  so  that  error  due  to  the  time 
discretization  would  be  negligible.  However,  in  the  present  algorithm 
the  stability  limits  permit  very  large  time  steps  so  that  the  time 
discretization  error  becomes  significant.  For  this  reason  it  would  be 
desirable  to  implement  a  second  order  accurate  version  of  the  new 
algorithm.  It  is  now  necessary  to  apply  the  new  method  to  modeling  the 
nonlinear  evolution  of  resistive  instabilities  in  realistic  tokamak  and 
reversed-f ield  pinch  equilibria. 
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FIGURE  CAPTIONS 
Fig.  1   Kinetic  energy  of  a  shear  Alfven  wave  due  to  a  three-dimensional 
perturbation,  with  w  -  0.^   and  At  -  0.1.   The  semi-implicit 
method  properly  simulates  the  wave  even  though  the  time  step  is 
four  times  the  usual  explicit  CFL  limit. 

Fig.  2  Kinetic  energy  of  the  wave  of  Fig.  1  except  with  At  =  0.5. 
This  time  step  is  more  than  twice  the  usual  shear  Alfven  CFL 
limit  for  this  20  mode  case.  The  frequency  is  correct  although 
damping  is  present  due  to  the  predictor-corrector  scheme. 

Fig.  3  Kinetic  energy  for  a  shear  Alfven  wave  produced  by  a  large 
amplitude  perturbation.  For  this  many  mode  case  the  time  step, 
At  =  0.1,  violates  the  usual  stability  conditions  not  only  for 
the  equilibrium  field,  B^,  but  also  for  the  large  amplitude 
perturbed  field,  B  .  Nevertheless,  using  the  semi-implicit 
method  stability  is  ensured  and  the  correct  frequency, 
to  =  1.0,  is  obtained. 

Fig.  4  Safety  factor,  q,  as  a  function  of  radius  for  a  cylindrical 
tokamak  equilibrium  which  is  unstable  to  an  m  =  2,  n  =  1 
resistive  instability. 

Fig.  5  Kinetic  energy  as  a  function  of  time  for  the  unstable 
equilibrium  of  Fig.  iJ.  The  time  step  exceeds  the  usual  CFL 
limit  for  the  shear  Alfven  modes  by  more  than  a  factor  of  ten, 
yet  the  growth  rate  is  accurate  to  within  a  few  percent. 
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Fig.  6  Radial  velocity  'profile  for  the  m  -  2,  n  -  1  mode  during  linear 
growth  from  the  equilibrium  of  Fig.   4. 

Fig.  7  J   profile  for  the  m  «  2,  n  -  1  mode  during  linear  growth  from 
the  equilibrium  of  Fig.   4. 

Fig.  8  Helical  flux  contours  for  the  case  of  Fig.   4  showing  the  m  =  2, 
n  =  1  island  near  the  saturation  of  linear  growth. 
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